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1.  Objective 


We  propose  to  develop  a  multi-scale  framework  to  simulate  and  predict  shear  initiated  reactions 
in  energetic  and  reactive  materials;  implement  this  technology  into  the  Combined  Hydro  and 
Radiation  Transport  Diffusion  (CTH)  hydrocode;  verify  the  algorithm;  and  conduct  preliminary 
studies  to  assess  this  novel  framework.  The  final  product  is  an  improved  predictive  capability 
implemented  into  the  CTH  hydrocode  developed  by  Sandia  National  Laboratories  to  allow  for 
the  prediction  and/or  simulation  of  high  explosive  (HE)  and  reactive  material  (RM)  reactions 
when  subjected  to  loads  that  result  in  shear  localizations. 


2.  Approach 


A  numerical  framework  for  nucleating  and  propagating  shear  localizations  already  exists  in 
CTH,  where  a  physical  criterion  for  the  initiation  of  shear  bands  has  been  implemented.  We 
propose  to  use  this  framework  to  track  shear  localizations  and  apply  the  outcome  of  our  multi¬ 
scale  approach  at  these  locations  to  simulate  the  shear  initiated  reactions  of  HEs  and  RMs.  A 
schematic  of  the  approach  developed  in  this  work  is  shown  in  figure  1.  It  consists  of  four  major 
components:  (a)  mesoscale  modeling,  (b)  converting  mesoscale  output  to  input  for  the  reactive 
bum  model,  (c)  linking  the  reactive  bum  model  output  to  continuum  level  simulation  of  shear 
initiated  reactions,  and  (d)  assessing  the  overall  approach  by  implementing  it  into  the  continuum 
mechanics  code  CTH. 
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As  depicted  in  figure  1(a),  each  mesoparticle  is  intended  to  represent  several  atomic-level  unit 
cells  of  the  crystalline  material.  The  dissipative  particle  dynamics  method  (1,  2)  provides  the 
mesoscale  input  to  the  reactive  bum  model.  Thermodynamic  properties  (pressure,  density, 
temperature,  and  energy)  from  the  mesoscale  simulations  are  generated  at  various  steady-state 
shear  rates,  providing  look-up  data  tables  for  use  in  the  CTH  modeling  (cf.  figure  1(b)),  such  that 
for  each  strain  observed  in  the  computational  domain  where  a  shear  localization  has  developed, 
the  state  variables  and  the  energy  spectrum  obtained  by  the  mesoscale  model  are  known  at  the 
continuum  level.  With  this  approach,  the  assumptions  made  in  the  continuum  mechanics  code 
with  respect  to  shear  localizations  are  circumvented.  Moreover,  by  incorporating  mesoscale 
level  calculations  into  the  continuum  mechanics  code,  more  accurate  estimates  for  the 
thermodynamic  state  within  the  localized  regions,  and  consequently  more  accurate  input  to  the 
reactive  bum  model,  are  obtained. 

A  sub-grid  model  is  developed  to  account  for  energy  dissipation  and  link  the  mesoscale 
temperature  estimates  to  the  continuum  level  in  a  relatively  cell-size  independent  manner  by 
calculating  a  volumetric  averaged  solution  of  the  extent-of-reaction  at  the  sub-grid  level  to 
represent  the  continuum  extent-of-reaction  (cf.  figure  1(c)).  This  extent-of-reaction  is  used  in  the 
equation-of-state  (EOS)  to  determine  the  amount  of  energy  release. 

The  sub-grid  calculations  are  performed  for  cells  that  are  identified  to  contain  shear  bands  at  the 
continuum  level.  A  numerical  framework  for  nucleating  and  propagating  shear  band 
localizations  already  exists  in  CTH  (5-5).  When  a  set  of  nucleation  criteria  is  satisfied,  a  shear 
band  is  formed  by  introducing  Lagrangian  particles  and  they  conform  to  local  planes  of 
maximum  shear  as  they  propagate  in  three-dimensional  space,  until  a  set  of  growth  criteria  is  no 
longer  satisfied  along  the  points  defining  its  boundary.  This  framework  is  used  to  track  shear 
localizations,  and  our  multi-scale  approach  is  applied  at  these  locations  to  simulate  the  shear 
initiation  of  energetic  and  reactive  materials  (cf.  figure  1(d)).  The  approach  is  implemented  into 
CTH  and  simulations  are  conducted  for  HMX  to  assess  the  algorithm.  Details  of  the  approach 
are  further  discussed  in  the  following  sections. 

2.1  Mesoscale  Methods 

The  energy-conserving  version  of  the  Dissipative  Particle  Dynamics  (DPD)  method  is  a  particle- 
based  mesoscale  method  that  conserves  both  momentum  and  energy  while  allowing  the 
mesoparticles  to  exchange  both  viscous  and  thermal  energy  (I,  2).  Further  details  of  the  DPD 
approach  can  be  found  in  the  original  papers  (I,  2)  and  the  report  for  our  Year-1  effort  (6). 

As  in  the  Year-1  effort,  steady  planar  shear  flow  was  induced  by  means  of  the  Lees-Edwards 
boundary  conditions  (7),  in  which  the  simulation  box  and  its  images  centered  at  (x,y)  =  (±L,  0), 
(±2 L,  0), ...,  are  taken  to  be  stationary,  while  boxes  in  the  layer  above,  (x,  y)  =  (0,  L),  (±L,  L ), 

(±2 L,  L ), ...,  are  moving  at  a  speed  mL  in  the  positive  x  direction,  where  mis  the  shear  rate. 
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Boxes  in  the  layer  below,  (x,  y )  =  (0,  -L),  (±L,  - L ),  (±2 L,  -L), ...,  move  at  a  speed  mL  in  the 
negative  x  direction.  A  system  under  such  conditions  is  subjected  to  a  uniform  steady  shear  in 
the  xy  plane. 

2.2  Mesoscale  Models 

Mesoscale  models  were  generated  for  both  HMX  and  a  nickel  (Ni)/aluminum  (Al)  system.  For 
the  HMX  case,  one  mesoparticle  was  chosen  to  represent  a  single  HMX  molecule  (figure  2), 
where  mesoparticles  interact  through  a  standard  pairwise  Morse  potential  (5).  Potential 
parameters  can  be  determined  by  analytical  solution  of  the  Morse  potential  expression,  where  the 
cohesive  energy,  density,  and  bulk  modulus  at  zero  temperature  and  pressure  are  inputted 
conditions.  The  inputted  cohesive  energy,  density,  and  bulk  modulus  can  be  taken  from  either 
experimental  or  higher  resolution  simulation  data.  Details  for  determining  the  Morse  parameters 
can  be  found  elsewhere  (P).  For  the  HMX  mesoscale  model,  good  agreement  was  found  for  both 
the  EOS  properties  (figure  3)  and  the  elastic  coefficients. 


Figure  2.  Mesoscale  model  representation  of  FIMX. 


a.D  D.1  0.2  D.3  0.4  D.5 

(V0-V)/Vtt 


Figure  3.  Pressure-volume  isotherms  (0  K)  for  an  atomistic  model  (black),  a 
theoretical  fit  (red),  and  our  mesoscale  model  (green)  of  FIMX. 
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For  the  Ni/Al  system,  one  mesoparticle  was  chosen  to  represent  four  face-centered-cubic  unit 
cells,  or  16  atoms,  where  mesoparticles  interact  through  a  Sutton-Chen  potential  (10).  The 
Sutton-Chen  potential  is  the  simplest  of  the  embedded-atom  type  potentials,  where  a  pairwise 
repulsive  interaction  is  combined  with  a  many-body  cohesive  term.  The  simple  analytical  form 
of  the  potential  has  contributed  to  its  popularity  for  simulating  metals  and  alloys.  Analogous  to 
the  determination  of  the  Morse  potential  parameters  for  HMX,  the  Sutton-Chen  potential 
parameters  can  be  determined  through  analytical  solution  of  the  potential,  where  the  cohesive 
energy,  density,  and  bulk  modulus  at  zero  temperature  and  pressure  are  inputted  conditions,  and 
again  these  property  values  can  be  taken  from  either  experimental  or  simulation  data.  Details  of 
the  Sutton-Chen  potential  and  the  approach  for  fitting  the  material  parameters  can  be  found 
elsewhere  (10). 

The  melting  behavior  of  these  models  must  be  accurately  captured  since  alloying  of  the  Ni/Al 
system  is  diffusion  driven.  Therefore,  further  fitting  of  the  model  parameters  to  match  the 
experimentally  determined  melting  temperature  was  performed.  For  both  the  Ni  and  A1  models, 
reasonable  agreement  was  found  for  the  EOS  properties,  as  shown  in  figure  4. 


V  [cmVmol] 


Figure  4.  EOS  property  comparison  between  atomistic  (MD)  and  mesoscale  (DPD)  models  for  Ni  and  Al. 

2.3  Sub-Grid  Models 

A  sub-grid  model  is  developed  to  account  for  the  dissipation  of  energy  and  allow  for  delayed 
reactions.  Initiation  phenomena  in  heterogeneous  energetic  materials  can  occur  when  the 
material  is  subjected  to  impulses  such  as  shock  waves.  Under  sufficiently  strong  shock 
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conditions,  these  shock  waves  can  evolve  into  self-sustaining  detonation  waves.  This  process, 
shock-to-detonation  transition  (SDT),  is  fairly  well  understood  at  the  phenomenological  level 
and  is  based  on  the  theory  of  hot  spot  formation  (such  as  void  collapse,  visco-plastic  heating, 
shear  band,  frictional  heating,  etc.).  However,  for  conditions  such  that  the  SDT  process  does  not 
occur,  initiation  of  the  HE  leading  to  explosion  can  take  place  if  the  material  is  sufficiently 
confined.  Such  analysis  was  shown  by  Frey  (9),  whose  study  indicated  that  both  pressure  and 
shear  rate  were  important  parameters  in  controlling  runaway  explosion.  In  another  study,  a  shear 
banding  mechanism  was  shown  to  provide  a  large  ignited  surface  area,  which  is  believed  to  be 
necessary  to  explain  shock  initiation  (10).  The  long  delay  between  the  time  at  impact  and 
explosion  can  typically  be  hundreds  of  milliseconds,  whereas  for  the  SDT  process  it  is  on  the 
order  of  microseconds.  Such  dominant  physical  features  suggest  energy  competing  processes  at 
the  micro-scale  level.  Therefore,  at  hot  spot  locations,  e.g.,  shear  surfaces,  the  possibility  of 
initiation  hinges  on  a  balance  (or  lack  thereof)  between  energy  producing  mechanisms  (visco¬ 
plastic  work,  shear  localization,  chemical  reaction,  etc.)  and  the  rate  at  which  the  energy  is 
transported  away  from  the  zone.  Energy  generated  from  shear  localization  within  the  narrow 
region  of  the  shear  band  width,  which  is  on  the  order  of  micrometers,  can  be  much  higher  than 
the  bulk  energy  in  the  typical  cell  size  of  continuum  simulations,  which  is  on  the  order  of 
millimeters.  The  large  temperature  gradient  within  a  continuum  cell  supports  the  need  to 
account  for  thermal  diffusion.  Figure  5  illustrates  temperature  profiles,  including  energy  release 
due  to  reaction,  as  time  progresses.  This  concept  forms  the  basis  of  our  sub-grid  model. 
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Figure  5.  Sheared  cell  and  its  temperature  profile. 
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Within  each  sheared  continuum  cell,  we  assume  that  the  shear  band  surface  is  located  at  the  mid¬ 
plane  (figure  5).  Such  a  geometric  assumption  leads  to  a  simple  and  tractable  solution  at  the  sub¬ 
grid  level.  Each  sheared  cell  is  subdivided  into  intervals,  and  time-dependent  equations  of 
temperature  and  species  are  described  as  follows: 


rST  lr 

pC—  =  k 


dt  d2x 


+  A Hm" 


pd-^  =  pDd-^-m" 
dt  d2x 


where  A  is  the  extent  of  reaction  at  the  sub-grid  level,  AH  is  the  heat  of  reaction,  D  is  the 

diffusion  coefficient,  k  is  the  thermal  conductivity,  and  the  volumetric  mass  production  rate, 
m"' ,  is  governed  by  an  Arrhenius-type  reaction: 


M-k) 


\  ~T/ 

.jexp  yT 


where  A  is  the  frequency  and  the  activation  temperature  is  Ta .  The  species  diffusion  term  is 

assumed  negligible  in  the  current  study.  Solutions  for  the  previous  equations  are  integrated 
using  an  implicit  method  with  the  initial  conditions: 

T(x,0)  =  Tmeso,  0  <x<wsb 

T{x,0)  =  TCe„,  Wsh  <x<U/2 

where  the  shear  band  temperature,  Tmeso ,  is  obtained  from  mesoscale  simulations  of  a  shear 
banding  process,  Tcell  is  the  continuum  cell  temperature,  and  wsb  is  the  shear  band  width. 
Boundary  conditions  are  given  as 


—  =0  tm&T{x  =  U/2)=Tcell. 


Within  each  sheared  cell,  a  volumetric  averaged  solution  of  the  extent-of-reaction  at  the  sub-grid 
level  is  transferred  to  the  continuum  level  as  the  continuum  extent-of-reaction.  This  extent-of- 
reaction  is  used  in  the  EOS  to  determine  the  amount  of  energy  release.  Such  an  approach 
ensures  consistency  in  our  internal  energy  calculation. 
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3.  Results 


Hydrostatic  simulation  was  carried  out  with  HMX  mesoscale  model  to  probe  the  range  of 
validity  of  our  model  as  well  as  to  verify  model  parameters.  Figure  6  shows  curve  fittings  of  the 
mesoscale  data  using  Hugoniot  relations.  Values  of  the  sound  speed,  the  Us-Up  slope,  the  bulk 
modulus,  and  its  derivative  are  within  ranges  reported  in  literature. 


Figure  6.  (a)  Hugoniot  fit  of  the  Us-Up  relation  and  (b)  Hugoniot  fit  of  the  P-V  relation. 


Analogous  to  our  Year-1  effort,  the  mesoscale  look-up  table  was  generated  for  a  range  of 
densities  and  shear  rates,  where  the  pressure  and  temperature  were  calculated.  As  an  example  of 
a  particular  data  point,  temporal  profiles  of  the  various  energy  contributions  for  the  HMX 
mesoscale  model  under  shear  is  shown  in  figure  7(a).  A  configuration  snapshot  of  the  HMX 
material  just  prior  to  the  release  of  elastic  strain  energy  is  given  in  figure  7(b),  where  evidence  of 
shear  localization  is  present. 


Figure  7.  (a)  Temporal  variations  of  various  energy  contributions  after  the  system  has  been  sheared  in  the 

mesoscale  simulation  and  (b)  a  configuration  snapshot  of  the  material  sample  just  prior  to  the  release 
of  elastic  strain  energy. 
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In  the  continuum  level  simulations,  a  block  of  material  is  subjected  to  pure  shear  by  imposing  a 
velocity  on  the  lower  half  of  the  block  to  verify  the  newly  implemented  approach  in  CTH,  as 
shown  in  figure  8(a).  Example  output  from  a  CTH  simulation  in  which  the  sub-grid  model  has 
been  implemented  is  shown  in  figures  8  (b)  and  (c). 


X,  Shearing 
surface 

(a) 


Extent  of  Reaction(2) 


40  - 
35  - 
30 

25  - 
20  - 
15  - 
1  0 
5 

0  - 
-5  - 
-10 
-15  - 
-20 


Time=2.53146x10-6 


-30 


(b) 


I  1.000x10* 
3.388*10-’ 
1.512x10-' 
S378*»-J 
2^85*»-J 
8386*10-* 
3455x10- 1 
1.343*10- 3 
5-223XKT'1 
2.031x10-* 
7397*10-* 
3.070*10-* 
1.194*10"* 
4.542x10-* 
1.805x10-“ 
7.017x10-' 
2.728*»-7 
1.061x10-' 
4.125x10"  8 
1.604x10-* 


(C) 


Figure  8.  Example  problem  to  verify  the  new  multi-scale  approach:  (a)  a  sliding  half  of  a  block  of  material,  (b)  the 
extent-of-reaction  at  2.5  ps,  and  (c)  the  extent-of-reaction  at  5  ps. 


4.  Conclusions 


We  developed  a  multi-scale  approach  for  simulating  shear  initiated  reactions  that  spans  from  the 
molecular  scale  to  the  mesoscale  (length  scale  of  material  heterogeneities)  to  the  full  continuum 
scale  (length  scale  of  the  weapon  system)  and  implemented  it  into  the  CTH  hydrocode  developed 
by  Sandia  National  Laboratory.  We  demonstrated  that  the  new  approach  allows  predictions  not 
previously  possible  for  energetic  materials  when  subjected  to  loads  that  result  in  shear 
localizations.  This  computational  tool  provides  a  novel  modeling  capability  that  opens 
previously  unavailable  avenues  by  bridging  the  gaps  between  multiple  scales  to  enable  improved 
predictions  toward  designing  armor  and  anti-armor  devices.  It  also  enables  the  development  of 
concepts  for  enhanced  survivability  and  lethality,  primarily  in  the  areas  of  insensitive  munitions, 
reactive  armor,  and  the  development  of  novel  concepts  and  designs  using  RMs. 

Accomplishments  in  this  year’s  work  included  the  development  and  implementation  of  a  simple 
and  efficient  approach  to  generating  various  mesoscale  models  for  HE  and  RM  systems.  We 
found  the  models  to  adequately  reproduce  the  EOS  properties  for  HMX  and  a  Ni/Al  system. 
Compared  to  Year  1,  the  approaches  used  to  develop  this  year’s  models  were  more  efficient 
while  still  maintaining  physically  relevant  parameters.  In  this  year’s  effort,  we  also  refined  our 


8 


sub-grid  model  from  the  previous  year.  Furthermore,  we  refined  the  multi-scale  links  between 
the  mesoscale  and  the  sub-grid  models  as  well  as  the  sub-grid  and  continuum  descriptions. 
Overall,  these  improvements  have  led  to  more  physically  pertinent  continuum  level  predictions. 


9 


5.  References 


1.  Bonet  Avalos,  J.;  Mackie,  A.  D.  J.  Chem.  Phys.  1999,  111,  5267-5276. 

2.  Mackie,  A.  D.;  Bonet  Avalos,  J.;  Navas,  V.  Phys.  Chem.  Chem.  Phys.  1999, 1,  2039. 

3.  Silling,  S.  A.  CTH Reference  Manual:  Three  dimensional  shear  band  model.  Unpublished. 

4.  Fermen-Coker,  M.  Numerical  Simulation  of  Adiabatic  Shear  Bands  in  H-6AI-4V Alloy  Due 
to  Fragment  Impact,  ARL-RP-91;  U.S.  Army  Research  Laboratory:  Aberdeen  Proving 
Ground,  MD,  March  2005. 

5.  Fermen-Coker,  M.  Implementation  of  Schoenfeld-Wright  Failure  Criterion  Into  a  3D 
Adiabatic  Shear  Band  Model  in  CTH;  ARL-TR  3284;  U.S.  Army  Research  Laboratory: 
Aberdeen  Proving  Ground,  MD,  September  2004. 

6.  Brennan,  J.  K.;  Tran,  L.;  Fermen-Coker,  M.  Physics  Based  Multi-Scale  Modeling  of  Shear 
Initiated  Reactions  in  Energetic  and  Reactive  Materials;  U.S.  Army  Research  Laboratory: 
ARL-TR  4753;  Aberdeen  Proving  Ground,  MD,  March  2009. 

7.  Lees,  A.W.;  Edwards,  S.  F.  J.  Phys.  C:  Solid  State  Phys.  1972,  5,  1921. 

8.  Allen,  M.  P.;  Tildesley,  D.  J.  Computer  Simulation  of  Liquids;  Oxford  University  Press: 
New  York,  1989. 

9.  Girifalco,  L.  A.;  V.  G.;  Weizer,  V.  Phys.  Rev.  1959, 114,  687. 

10.  Sutton,  A.  P.;  Chen,  J.;  V.  Phil.  Mag.  Lett.  1989,  61,  139. 

11.  Frey,  R.  Seventh  Symposium  (International)  on  Detonation,  Naval  Surface  Weapons  Center, 
NSWC  MP  82-334.  White  Oak,  MD,  1981. 

12.  Howe,  P.;  Frey,  R.;  Taylor,  B.;  Boyle,  V.  Sixth  Symposium  (International)  on  Detonation, 
Office  of  Naval  Research,  ACR-21 1,  Arlington,  VA,  1970. 


10 


6.  Transitions 


Transitioning  this  modeling  capability  to  the  FY1 1  program  goals  of  the  Institute  for  Multi-Scale 
Reactive  Modeling  of  Insensitive  Munitions  (MSRM)  is  planned. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


A1 

aluminum 

CTH 

Combined  Hydro  and  Radiation  Transport  Diffusion 

DPD 

Dissipative  Particle  Dynamics 

EOS 

equation-of-state  properties 

HE 

high  explosive 

MD 

molecular  dynamics 

MSRM 

Multi-Scale  Reactive  Modeling  of  Insensitive  Munitions 

Ni 

nickel 

RM 

reactive  material 

SDT 

shock-to-detonation  transition 
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